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Abstract:  We  prove  that  for  solutions  to  the  two  and  three  dimensional  incompress¬ 
ible  Navier-Stokes  equations  the  minimum  scale  is  inversely  proportional  to  the  square  root 
of  the  Reynolds  number  based  on  the  kinematic  viscosity  and  the  maximum  of  the  veloc¬ 
ity  gradients.  The  bounds  on  the  velocity  gradients  can  be  obtained  for  two  dimensional 
flows,  but  have  to  be  assumed  in  three  dimensions.  Numerical  results  in  two  dimensions 
are  given  which  illustrate  and  substantiate  the  features  of  the  proof.  Implications  of  the 
minimum  scale  result  to  the  decay  rate  of  the  energy  spectrum  are  discussed. 
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1  Introduction 

We  consider  solutions  to  the  incompressible  Navier-Stokes  equations 
ut  +  uua  +  vu„  4-  wnz  +  Vp  =  mAu,  v  >  0, 

Um  +  vt  +  wt  =  0, 

on  a  2ir-periodic  square,  where  u  =  (u,v,w)  is  the  velocity  vector,  p  the  pressure 
and  i/  the  kinematic  viscosity.  Solutions  of  the  Navier-Stokes  equations  for  small 
viscosity  are  usually  turbulent;  such  flows  posses  a  lot  of  structure  in  both  space 
and  time.  The  viscosity  of  the  fluid  controls  the  level  of  turbulence  within  a  flow 
by  affecting  the  energy  dissipation.  As  the  viscosity  is  decreased  the  size  of  the 
smallest  features,  or  scales,  diminishes.  The  relation  between  the  viscosity,  the 
minimum  scale  and  the  total  energy  dissipation  is  of  fundamental  interest  for  the 
understanding  of  turbulence. 

The  mathematical  theory  for  the  Navier-Stokes  equations  is  not  complete  for 
three  dimensional  flows:  the  global  regularity  is  not  known  and  no  global  bound  for 
the  velocity  gradients  is  available.  However,  both  these  results  are  known  in  two 
space  dimensions.  Whether  the  results  from  two-dimensional  flows  are  of  physical 
relevance  is  open  to  discussion  since,  in  the  absence  of  viscosity,  flows  in  two  di¬ 
mension  conserve  both  energy  and  enstrophy,  while  in  three  dimensions  only  the 
energy  is  conserved.  Nevertheless,  results  on  two-dimensional  turbulence  may  be  of 
significance  for  large  scale  oceanographic  and  atmospheric  motions. 

Assuming  global  regularity,  we  relate  the  minimum  scale  of  the  flow  to 
the  global  bound  of  the  velocity  gradients.  Our  main  result,  precisely  stated  in 
theorem  (2.1),  is  that  the  minimum  scale  is  essentially  no  smaller  than 

By  comparison,  a  commonly  accepted  minimum  scale  for  two  dimensional  flows, 
Ajx) i  (see  Lilly  [19],  Orssag  [20]),  is  based  on  the  total  dissipation  rate  of  the  en¬ 
strophy  per  unit  volume.  The  enstrophy  is  defined  as  the  square  of  the  Lj-norm  of 
the  vorticity.  From  dimensional  arguments  it  follows  that 


where 


1  =  21 >f  JwU'  +  WU'd'dy 

is  the  total  rate  of  enstrophy  dissipation  per  unit  volume  and  £  is  the  vorticity. 

In  three  space  dimensions  the  corresponding  minimum  scale  is  the  Kolmogoroff 
dissipation  scale  [15] 

where 

«  =  2 "J  J  J  IK  II*  +  IK II*  +  IK  II 'dzdyiz, 

is  the  total  rate  of  energy  dissipation  per  unit  volume. 

The  estimates  for  the  minimum  scale  can  be  used  to  determine  the  decay  rate 
of  the  energy  spectrum,  assuming  that  a  power  law  does  in  fact  exist. 

From  our  results  in  two  dimensions  we  conclude  that  the  energy  spectrum, 
E(k),  behaves  like  k~a  when  there  is  a  maximum  rate  of  enstrophy  dissipation  in 
the  flow.  The  k~a  power  law  is  in  accordance  with  the  Batchelor-Kraichnan  theory 
of  enstrophy  cascade  [3]  [16].  The  high  rate  of  dissipation  can  not  remain  for  long 
times  without  the  flow  disappearing.  Indeed,  numerical  experiments  show  that  the 
solutions  rearrange  themselves  into  organised  structures  which  dissipate  enstrophy 
at  a  much  smaller  rate.  Saffman’s  work  [22],  which  predicts  a  power  law  Jfe~4, 
seems  to  describe  the  behavior  of  the  system  at  this  later  stage  of  evolution.  Our 
theory  does  not  predict  the  power  law  but  only  relates  it  to  the  rate  of  enstrophy 
dissipation;  the  k~ 4  law  would  correspond  to  rj  of  order 

In  three  space  dimensions  there  is  no  a  priori  bound  for  IDul^.  However,  if 
we  assume  that 

\Dn\oo 

then  when  the  energy  dissipation  rate  e  is  of  order  one,  we  obtain  the  Kolmogoroff 
power  law,  E(k)  =  fc-*/8,  and  the  Kolmogoroff  scale  Ami„  =  Asp  =  i/8^4. 

Some  of  the  first  calculations  on  two-dimensional  turbulence  were  performed  by 
Lilly  [19],  Fox  and  Orszag  [9],  Herring,  Orszag,  Kraichnan  and  Fox  [12],  Fornberg  [8], 
and  Barker  [1],  among  others.  More  recent  computations  on  meshes  of  up  to  1024  x 
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1024  points  are  described  in,  for  example,  Brachet,  Sulem  [7],  Brachet,  Meneguzzi, 
and  Sulem  [6],  Herring,  and  McWilliams  [11]  and  Benzi  et  al  [4].  In  some  cases  the 
small  viscosity  limit  of  the  equations  was  approximated  by  the  continuous  removal 
of  the  high  frequency  Fourier  coefficients  [8].  In  other  cases  the  true  dissipation 
term  is  integrated  although  some  extra  smoothing  of  high  frequencies  is  sometimes 
required  to  suppress  the  growth  of  aliasing  errors  [l] .  Another  approach  is  to  replace 
the  viscosity  term  by  a  super-vitcosity,  that  is  a  higher  power  of  the  Laplacian 
operator  [7],  [11].  This  operator  allows  simulations  with  a  formal  viscosity  which  is 
much  smaller.  The  minimum  scale  is  nevertheless  comparable  to  the  computations 
presented  in  this  paper. 

The  numerical  simulation  of  three  dimensional  flows  is  still  limited  by  the 
power  of  current  computing  machines.  Currently,  the  largest  three  dimensional 
simulations  seem  to  have  been  performed  on  128s  meshes.  However,  by  exploiting 
the  symmetries  of  the  Taylor-Green  problem,  Brachet  et  al.  were  able  to  effectively 
solve  with  a  256s  resolution  [5].  They  find  the  slope  of  the  spectrum  to  be  least 
steep  when  the  rate  of  energy  dissipation  reaches  a  maximum.  The  numerical  results 
seem  to  agree  at  this  point  with  the  Kolmogoroff  scale.  For  further  references  on 
three  dimensional  computations  see  the  review  article  by  Hussaini  and  Zang  [14], 

We  restrict  ourselves  to  two  dimensional  simulations.  Our  numerical  approach 
has  been  to  attempt  to  faithfully  solve  the  viscous  Navier-Stokes  equations.  The 
computations  were  performed  using  the  pseudo-spectral  method,  Kreiss  and  Oliger 
[18],  and  Orszag  [21].  There  is  no  extra  viscosity  added  to  the  numerical  simulation 
through  smoothing  or  chopping  of  the  high  frequencies,  although  the  fourth-order 
predictor  corrector  time  integrator  produces  a  small  amount  of  it.  Numerical  simu¬ 
lations  are  used  to  confirm  the  theoretical  estimates  and  to  show  that  the  estimates 
can  be  achieved  for  certain  initial  conditions.  Results  are  shown  for  the  time  de¬ 
velopment  of  a  flow  which  initially  is  maximal  dissipative.  We  also  show  results  of 
forced  problems.  In  this  case  there  is  no  easy  a  priori  bound  for  the  maximum  norm 
of  the  vorticity.  We  found  numerically  that  the  forcing  should  be  proportional  to 
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the  viscosity  in  order  to  obtain  order  one  velocities.  More  numerical  work  on  this 
subject  is  still  necessary. 

In  section  2  we  present  the  analytical  results.  In  section  3  we  present  numerical 
computations  in  two  space  dimensions  which  substantiate  and  illustrate  various 
features  of  the  proof.  Finally,  in  section  4  we  discuss  the  implications  of  the 
minimum  scale  result  to  the  decay  rate  of  the  energy  spectrum. 

2  Analytical  Results 

In  this  section  we  will  prove  some  results  about  the  rate  of  decay  of  the  Fourier 
coefficients  for  solutions  of  the  incompressible  Navier-Stokes  equations, 


Uj  +  unB  4-  vuy  +  wut  +  Vp=  i/Au,  v  >  0, 
UB  +  Vy  +  Wt  =  0, 


(2.1.) 

(2.1b) 


on  the  region  12  =:  {0  <  x,y,z  <  2ir}  and  for  t  >  0.  We  assume  that  u  = 
(u(x,  t),  t;(x,  t),  w(x,  t ))  is  2ir-periodic  in  x  =  (*,  y,  z). 

At  t  —  0  we  give  the  initial  data 

u(x,  0)  =  u0(x)  ,  V  •  u0  =  0. 


For  simplicity  we  assume  that 


u0dx  =  0, 


which  implies 


I  u(x,  t)dx  =  0  ,  for  t  >  0. 

J  n 


We  assume  that  (2.1)  has  a  bounded  solution  for  all  times  and  want  to  show 
that  the  smallest  scale  is  essentially  proportional  to  (»//|  £>111^  )x  .  Here 


=  suplPuU  and  \Du\n  =  sup(|0u/0*|,  |0u/0y|,  |0u/0z|) . 

I  X 


In  general  let 


Dpu  = 


dxp'  dyp^dz**  ’ 


mm 
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denote  any  derivative  of  u  of  order  p,  where  p  =  px  -f  p2  +  p3. 

We  are  interested  in  the  case  when  v  <  1  and  >  const  >  0.  Let  us 

further  assume  that  for  every  natural  number  p  there  is  a  constant  Cp  such  that 
the  initial  conditions  satisfy  the  bounds 

I^Too 

U P 

(2.3) 

Here  |f|  oo  —  supxefj  |f(x)|  denotes  the  maximum  norm  and 

(f,g)=  [t  gdx  ,  ||f||2  =  (f,f)  , 

J  n 


max  Hf(0)  =:  max  ( 
o<,< r  ,v  7  o<j<f ' 


lf^(.o)HI)  +  |||^(..o,|p, 


the  1/2-scalar  product  and  norm.  Then  we  prove 

Theorem  (2.1):  We  assume  (2-3)  holds  and  develop  u  into  a  Fourier  series 
«(*.<)  =  (M)e*'<k,x>,  k  =  {kuk3,k5). 

k 

For  every  natural  number  j  and  any  real  number  a  >  0  there  are  constants  K  and 
K  which  depend  on  j,  a  and  Ci  (l  =  l(j,a)),  such  that 


and 


sup|u(M)|l 

«>o 


<  K 


r— - — -j  +  l  +  a 

IPiiL 

J/i-i+Q  |fc|2j 


sup(u(fc,<)|2  <  K 
«>  o 


|Ou| 


;+a 

oo 


l/J+a  \k\2>' 


The  estimate  of  the  theorem  can  be  rewritten  in  the  form 


sup  ju(fc,t)|2  <  K 

l>0 


ij 


We  see  that  the  spectrum  becomes  vanishingly  small  once  |fc|  (|L>*1|00/*/)1^J  with 
|u(M)l  decaying  faster  than  any  power  of  (|£>u|00/t/)1/2|At|_1.  It  is  natural  to  define 
the  minimum  scale  of  the  flow  to  be  proportional  to  (*//|-Du|0C)1/2. 


In  three  space  dimensions  there  is  no  a  priori  bound  for  l-Du^.  One  can 
speculate  what  the  right  order  of  magnitude  is.  For  example,  if  we  assume  that 

|0«L  ~ 

we  obtain  _ 

(>0  u 

which  corresponds  to  the  Kolmogoroff  scale  [19]  of  Am,-n  =  i/s/4. 

In  contrast,  for  two  space  dimensions  an  a  priori  bound  for  IDu^  can  be 
obtained.  The  vorticity,  £  =  ut  —  vm,  which  satisfies 

&  4-  -f  v£y  =  vA£, 

obeys  the  maximum  principle 

Moo  =  sup  |£(x,  f)|  <  sup  |£(x,  0)|. 

x,<  x 

Therefore,  assuming  that  the  initial  data  satisfy 

sup  |£(x,  0)1  <  1,  (2.4) 

X 

we  shall  prove  that  [Uu^  is  essentially  bounded  independent  of  v,  in  the  sense 
that  for  every  /3  >  0  there  is  a  constant  K\  —  Ki(/3)  such  that 

|0u|oo  <  Kxv~*. 

Thus  if  the  initial  data  have  derivatives  of  order  one  then  the  smallest  scale  is 
essentially  of  the  order  vx^ . 

Our  proof  is  also  valid  for  Burger’s  equation.  In  this  case  one  can  prove  that 
\Du\  oo  ^  const  u  sup  |  Z?u(* ,  0)| oo • 

a 

Thus  our  result  predicts  that  the  minimum  scale  is  of  order  i/”1.  This  bound  can 
be  attained  in  the  presence  of  shocks. 


2.1  Estimates  for  p  <  3 

From  now  on  we  shall  assume  that  an  estimate  for  |  Z?u|  ^  exists.  Integration 
by  parts  give  us  the  basic  energy  estimate 


where 


Since  by  assumption  (2.3),  ||u(-,0)||2  <  const  ,  it  follows 

2v  f  H\(t)dt  <  ||u(-,0)||2  <  const  and  J|u(-,t)j|2  <  ||u(*,0)j|2  <  const  .  (2.6) 
Jo 

Now  differentiate  (2.1).  For  any  first  space  derivative  Du  we  obtain 

+  I  =  -u(l|Ou.||’  +  IIBn.ll1  +  IIOu.ll1), 

where 

I  =  (Du,  D(uua  +  vut  -f  toUj  ))  =  //+  III, 

II  =  (Du,  uDua  +  vDuv  +  wDuz), 

III  =  (Du,  Du  u9  +  Dv  Uy  +  Dw  u«). 

Integration  by  parts  and  V  •  u  =  0  shows  that  II  =  0.  Again  by  integration  by 
parts  we  obtain 

III  <  const  IJDuloo-ff?. 


Therefore 


that  is 


2  fit  II^UH2  - con,t  \Du\°°H*  ~  +  llI>u»llJ  +  llr>u*ll2). 


const  \Du\nH\-vHl 


Integrating  the  last  inequality  with  respect  to  t  gives  us, 

H*(t)  <  Hl(0)  +  const  IDul^  f  J/2(r)dr  -  2v  f  H\(t)<1t. 

Jo  Jo 


.  i’.  i'i  *'k  ixi  i  i  •«*  #-«*.♦  ti* 


Therefore  by  (2.3)  and  (2.6) 


H\it)  <  const  -^-■qg-  and  v  [  H2(t)dt  <  const  ■ — — . 

v  Jo  v 


For  the  second  derivatives  we  obtain 


where 


^||Z>*u||2  4  /  =  -«/(||I>2uB||2  +  ||Z?2u,i|2  4  IIUVII2), 


I  =(D2u,  D2(uu*  4  4  wnt))  =  II  4  III  4  IV , 

II  =(D2u,uD2Ug  4  vD2  u„  +  iu£>2ut)  =  0, 

III  =2(D2u,  DuDua  4-  DvDuv  4  DwDut)  <  const  |.Du|0O.H'2, 

IV  ={D2 u,  D2u  u,  4  D2v  u#  4  D2w  uz)  <  const  \Du\ooHl. 


Therefore, 


2qIH *  -  conat  \Du\°°Hl  ~  "Hi 


Integrating  the  last  inequality  with  respect  to  i  and  using  (2.3),  (2.7)  gives  us 

H\[t)  <  const  ^UJ°°  and  v  [  H%{t)dt  <  const  ^  ( 

v2  J0  v2 


For  the  third  derivatives  we  obtain 

i(Z>V  u),  +  I  =  -r(||Bsu.ll!  +  l|Bsu,ll!  +  II 0s". II2). 

where  by  Leibniz’s  rule 

/  =  (D9n,  D9( uu„  +  t)Uj  +  u>ur))  =  II  4  III  4  IV  4  V, 

II  —  (D3u,uD9um  4  vD9 u,  4  wD9ut)  =  0, 

III  =  3  (Dsu,Z?u  D2um  4  Dv  D2Uj  4  Dw  D2ut)  <  const  \Du\ooHl , 

IV  =  3  (D9u,D2u  Dus  4  D2v  Dut  4  D2w  Dut)  <  const  {DulooH+Hi, 
V  =  (D9 u,  D9 u  u,  4  D9v  uv  4  D9w  ut)  <  const  IDul*,#/. 


'  A  vS  •-> -AS  V  *>*/v\rV  .Sv.-  v.  /'.  /  . 
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Therefore 


-~Hl  <  const  \Du\niH}  +  H<H2)  -  uH\ , 

2  at 

H  2  j 

<  const  \Du\ooHl  +  const  [Duj^—  - 

1/  JL 


thus  using  (2.3),  (2.7)  and  (2.8) 


_ s  _ s 

Hl(t)  <  const  - — j-55-  and  v  J  H\(t)dt  <  const  - — J22-. 


2.2  The  estimates  for  general  p 

We  now  prove  theorem  (2.1)  for  arbitrary  p.  First  we  obtain  energy  estimates 
for  HP  in  terms  of  1  <  j  <  max(l,[(p  -  l)/3]).  These  estimates  are  used 

to  obtain  bounds  for  |.D,'u|00  in  terms  of  iDuj^.  Finally  improved  estimates  are 
obtained  for  Hp ,  l^ul^  and  | JO-*  u| ^  using  interpolation  inequalities;  the  theorem 
then  follows.  We  start  with 

Lemma  (2.1):  For  every  p  there  is  a  constant  Kp  such  that 
|(Z?,’u,U,’(uuB  +  vv^  +  urn*)  | 

/  m«*(l,[(|>-l)/3])  \ 

<  Kp  I  \D*UHl  +  Hp +1  J2  l^uU Hp_jj  . 


[i]  =  largest  integer  <  x  and  )DJu|oo  =  max  |Dfcu|. 


Proof:  We  need  to  estimate  expressions  of  the  form 

{Dpu,D*-ku  Dkum  + D*-kv  D'up  +  Dv~h-w  Dhnt)  for  A  =  0,...,p-1. 

We  integrate  by  parts  to  decrease  the  order  of  Dv u.  In  doing  this  the  order  of 
Dr~hu,  Dr~hv  and  Dr~hw  or  the  order  of  Dk uB,  Dh uv  and  Dkut  will  increase. 
For  each  new  term  generated  through  integration  by  parts, 

(D’u,  D^u  Dqium)  +  (I^u,  D^v  D" uF)  +  (£>’ u,  D"w  D" u,), 


we  can  decrease  the  order  g  and  increase  one  of  q2  or  gs  until  one  of  the  following 
conditions  is  satisfied: 

(1)  q  -  1  <  gs  <  q  and  q2  <  q 

(2)  q  <  qj  <  q  4  1  and  gs  <  g  —  1 
Note  that  in  case  (1)  if  gs  =  g  then 


{Dq u,  Dq,uDqum  4  Dq,vDquv  4  Dq,wDqnx)  =  0. 


It  follows  that 


I  =  (Dv u,  Dp(uum  +  vUy  4-  tuu,)  ) 


can  be  written  as  a  sum  of  terms 

A  :=  (Dq u,  D,p~*q+l u.  P9_1ua  +  D2p'wv  Dqlut  4  Dlr~,q+1w  Dq~l uz), 

1  <  2p  -  2g  4  1  <  g, 

B  ■-  {Dqu,  Dq+iu  D2p~2q  1  ua  4  Dq+1v  D2p-2q~1uv  4  Dq+1w  D2p~2q~1uz), 

1  <2p-2g-l  <g-l, 

C  :=  (Dqu,  Dqu  D2p~2gux  4  Dqv  D2p~2guv  4  Dq w  D2p~2quz), 

0  <  2p  -  2g  <  g  -  2. 

Expression  A:  If  2p  —  2g  4  1  =  1  then  the  estimate  follows,  otherwise  integration 
by  parts  is  applied  to  expression  A, 

A  =  {Dp-qu,Dp~q+l  ( DquDq~1nx )  )  4  {Dp~qv,  Dp~q+l  {DquDq-1uy)  ) 

4  ( Dp-qw ,  Dp~q+1  ( DquDq~1ux )  ) 

to  reduce  the  order  of  the  factors  D2p~2q+1u,  D*p~2q+lv  and  D2p~2q+lw.  In  this 
way  we  can  write  A  as  a  sum  of  terms 

(Dq'u1D*-qu  Dq,~luw  4  Dp~qv  Dq’~lur  4  Dp~qw  Dqi~lnz), 


where 


p  —  g  >  1,  <  P  4  1,  and  gi4gj=p4g4l. 


r-  _Nl.>  .*»  A  .V  \  a*  V  A  A  %tV>.  %,VA 


Also,  l<2p-2g+l<g  impbes 


p  —  q  <  max 


(i.li(P-i)])- 


The  required  estimates  can  be  obtained  in  the  following  cases: 

(i)  If  p  -  q  =  1  then  either  gi  =  qt  =  p  or  one  of  the  q j  =  p  +  1  and  the  other  is 
equal  to  p  —  1 . 

(ii)  One  of  the  q,  is  equal  to  p  +  1. 

When  neither  (i)  nor  (ii)  is  satisfied  then  g;-  <  p  +  1  and  p  —  q  >  1,  and  we  can 
reduce  Dp~q  further.  This  shows  that  A  can  be  estimated  in  the  desired  way. 
Expression  Bi  Correspondingly,  by  reducing  the  order  of  D2p~2q~lnm,  D2p~2q~1uy 
and  D2p~2q~l ut,  B  can  be  written  as  sum  of  terms 

( Dq'u,Dq'u  Dp-q-'um  +  Dqiv  Dp-q~luv  +  Dqiw  Dp~q~'  uz)  :=  II, 


where 


9j  <  P  +  1  >  9i  +  92  =  P  +  9  +  1  and  p  -  q  -  1  >0. 


Also  2p  —  2q  <  q  implies  0  <  p  -  q  -  1  <  [|(p  -  3)].  II  can  be  estimated  in  the 
following  two  cases: 

(i)  If  p  -  q  -  1  =  0  then  9i  =  p  +  1,  92  =  p  -  1  or  qx  —  q2  =  p,  and 

II  <  conat  iDuloo Hp+iHr-i,  or 
II  <  const  iDulooH*. 

(ii)  If  91  =  p  +  1  or  g2  =  p  +  1. 

Otherwise  q,  <  p  +  1  and  p-q  -  1  >0  and  hence  we  can  diminish  p  -  q  -  1  further. 
Therefore  we  obtain  the  desired  estimates  for  B. 

Expression  C:  Integration  by  parts  allows  us  to  diminish  the  order  of  D2p~2q u 
obtaining  terms  of  the  form 

(Dqi u,  D"3u  Dp-q-'um  +  Dqiv  Dp-q-'ut  +  Dq’w  Dp~q~1ut) 


where 


qj  <  p  +  1  ,  9i  +  92  =  p  +  9  +  1  and  p  -  9  -  1  >  0. 


Using  the  same  argument  as  for  B  we  obtain  the  desired  estimate.  This  proves  the 
lemma. 


Differentiating  (2.1)  p-times  gives  us 

(Dpu,  Dp u)  +  (Dpu,Dp(uum  +  uu,  +  urn,)  ) 

2  at 

- -►(lie'll. !|,+l|B'n,H,  +  l|B'n, ||’) 

Therefore  we  obtain  from  lemma  (2.1) 

/  maxCMfr-1)/8])  \ 

—  H*  <const  f  \Du\ooH*  +  Hp+i  |^u|oo^fp-;  I  —  2uHp+l 

(  m**(l,[(p-l)/S])  \ 

|Du|oo^*  +  -  ^  \D’u\loHp-j  I  ~  uHp+i 


Using  the  notation 


the  last  inequality  implies 


Lr  =  I"  H*(t)dt 
Jo 


(rjr~ r»  .  m«x(l,[(p-t)/Sj) _  \ 

H’(t)  <  conit  l!£tt+|5=[„I,+  I  •£  iU>uLi,-i)  (2  10) 


Lp+i  <  const 


TBSiL  ^  ieu| 


m«x(l,[(p  — 1)/3])  \ 

+  t  E  [1") 


To  begin  with  let  us  obtain  estimates  for  Hp  and  Lp+ 1,  for  p  —  4, . . . ,  7.  In  most 
applications  this  is  all  what  is  needed.  From  our  previous  results  we  know  that 

|DujS 

Ls  <  const  - — ^  and  T4  <  const  - 7 — . 


Therefore  (2.10)  and  (2.11)  give  us 


H\{t)  <const 


Li  <const 


+  leoLt.  +  i|o»lLii) 
+  !^Wi4  + 

\  I/S  V  v2 


D11 

<  const  - — j-92-, 


1 


<  const 


The  estimates  for  H*,  I8.  Hi  and  Lj  follow  in  the  same  way. 
For  p  =  7  we  obtain 


Hj{t)  <  const 


(¥■ 


+  {Dul^L-r  +  i|I?u|^Te  +  i|D2u|^Is  j  • 


Thus  we  have  to  estimate  iD’ul^.  By  Sobolev  inequalities 

10’ulL  <  rf,’+s  +  ~-Hl 

and  we  obtain  for  e  =  Hp/Hp+ s 


\JPn\l,  <  const  Hp+iHp. 


(2.12) 


In  particular 


|.D8u||0  <  const  H^Ha  <  const 


Now  we  use  the  usual  interpolation  inequality  (see  for  example  [10]) 

|X>Ju^  <  e\Dvu\l0  +  const  1  <  j  <P 

which  gives  us  for  e  =  (|Du|*0/|D,>u|*0)^,,''','^,>_1J 

|U'u|^  <  const  (2.13) 


For  j  =  2,  p  =  3  we  obtain 


|D!u|L  <  “mil  (|0*u|5.)1/,(|Cn|J0),,J  <  const 


(2.14) 


Therefore 

r— - j-7+1/4  —-7+1/4 

H](t)  <  const  and  I,  <  const  ^  —  • 

Using  (2.14)  we  can  in  the  same  way  estimate  H\  and  Hg.  These  estimates  are 
not  as  sharp  as  required  by  theorem  (2.1).  To  obtain  the  required  estimates  we 
have  to  estimate  /f  2  for  general  p  and  then  improve  the  lower  order  estimates  using 
interpolation  inequalities. 


Using  (2.10)  and  (2.11)  recursively  and  the  estimates  for  Hp  and  Lp +1  for  p  <  6 


we  obtain 
Lemma  (2.2): 


h;+1  < con.,  Y, 


(2.15) 


ii-ik 


where  the  sum  is  over  all  {k,  ji,p,,qi}  satisfying  the  constraints 

h  h 

Tk  =  p  +  1  -  ~  1)  =  2fc  +  g„ 

i=l  i=0 

1  <  Ji  <  [Pi/3]  ,  Pi+i  =pi  ~  ji  ~  1  ~qt  ,  0  <  g,  <  pi  -  ji  for  i  =  l,...,fc 

0  <  9o  <  p  +  1  ,  qu  =  pu  ~  jh  ,  pi  =  p  -  go  ,  fc  >  0. 

Note  that  ti,  =  p  +  1  for  p  <  5  . 

Proof:  We  first  obtain  estiamtes  for  Lp+2  from  (2.11).  H*+l  is  bounded  by  v  times 
this  estimate  for  Lp+ j. 

_  -  _  m»x(l,[p/S]) 

L,+,  <  const  (  /if*'  +  (I^T* /v)L,+ 1  +  £  |L»u|^L,,+  ,_;  ! 

<  const  (  A  +  B  +  ,_1  ] 

c 

where  we  have  labeled  the  three  terms  on  the  right  hand  side  as  A,  B  and  C.  In  the 
recursive  reduction  of  Lp+2  we  must  consider  all  possible  terms  which  may  arise; 
at  each  new  stage  one  must  consider  the  effect  of  using  expression  A,  B  or  C.  In 
the  general  case  one  chooses  B  go-times  followed  by  term  C  with  j  =  ji ,  then  B 
gi-times,  C  with  j  =  j2  and  so  on  until  finally  finishing  with  term  A  or  B  and  using 
the  estimates  for  Lp  with  p  <  6.  In  this  fashion  the  general  term  will  be 

\DuC,  \Dh*\L 

|/90  g/J  |/fll  j/1  +  l 


r*  n; 

'<»  ,.1^  +  lk  +  l 


We  define  r*  =  2k  +  £g<.  The  constraints  on  ji,  qi  and  p<  follow  from  the  manner 
in  which  the  general  term  was  obtained. 


Now  we  use  (2.13)  to  estimate  |i?'Mu|(X1  in  terms  of  and  |Z?u|  , 


<  const  Y.  ^ 

Jl  -Jk 

,  ri_s>  rrr-rl+Ptri-S)/^-*) 

<  const  £15=151  - - 


r*-2A 


Therefore  by  (2.2)  and  Sobolev  inequalities 

<  con.l  (HU,+«Ul) 

<  con»t  H*+  j 

<  const  max  !^! 


l+J>(r»-3)/(j>-2) 


— - - - rl  I  |DU 

Z?,’“1u!  <  const  max  I  — 


It  follows  that 


^  ^  J  • 

One  more  application  of  the  interpolation  inequalities  gives  for  1  <  j  <  p  —  1 

r7T-I> 

IZ^L  <  conit  m«  j  ■ 

In  order  to  obtain  our  final  estimate  we  need  to  show  that  that  t  =  min*  t* 
tends  to  infinity  as  p  tends  to  infinity.  Recall  that 

k  h 

n,  =  p  +  1  -  X)(;<  -  i)  =  2fc  +  ]Tgi, 

i=l  i= 0 

1  <  ji  <  [pi  /  3]  ,  Pi+i  =Pi  -  ji  -  1  -9i  ,  0  <  qi  <  pi  —  ji  for  t  =  l,...,fc 


o  <  9o  <  P  +  1  »  9fc  =  Ph  -  jk  ,  Pi  =  p  -  90  .  fc  >  o. 


From  ji  <  [pi/3]  <  Pi/3  it  follows  that 


This  last  expression  can  be  written  as 


*  >  1  + 


<.(^1 


Since  g*  =  Ph  -  jh  >  (2/3 )ph  we  have 


^  2  ■ _ ■ 

=  2k  4-  y] qi  >  2k  +  -ph  +  qj. 


Now  consider  the  two  cases 

0)  e?;„i(V2),»  >  (p + *)/* 

(a)  £&(»/»)'«  <  (p  +  3)/2 

In  case  (i)  it  follows  that 


£(3/2)fc_19.  >  £(3/2)*' gt-  >  (p  +  3)/2 

<=o  i=0 

fc-i 

=*  X>  >(2/3)fc~1(p  +  3)/2 

i=0 

=>  t*  >  2*  +  (2/3)*~I(p  +  3)/2 
Minimizing  this  last  expression  with  respect  to  k  (>  0)  gives 


T),  >  const  log(p)  for  some  const  >  0. 


In  case  (ii)  we  have 


‘-1+logs"(2(STj)) 

=>  ii  >  2  +  21ogs/1(^-ti_)  +  (2/3)pfc 
Since  0  <  p*  <  p  +  1  it  follows  that  as  a  function  of  p*,  the  above  expression  has 


the  bound 


Tfc  >  const  log(p)  for  some  const  >  0. 


Hence  r  =  min*  >  61og(p)  as  p  tends  to  infinity,  for  some  constant  6  >  0.  Thus 


.>  V  V 


Theorem  (2.2):  For  every  j  >  1  and  any  a  >  0,  we  can  choose  p  tufficiently  large 


so  that 


j  D*  u  |  <  const 


|Z>u 


ri  +  l+a 


(2.16) 


v  j-l+a  ’ 

where  the  constant  depends  on  p  and  Cr+ 1  introduced  in  the  estimates  for  the  initial 
conditions  (2.3)  . 

By  using  (2.16)  in  (2.15)  we  obtain 

Theorem  (2.3):  For  every  j  >  1  and  any  a  >  0,  we  can  choose  p  sufficiently  large 
so  that  + 

ff’(0  <  const  ,  (2.17) 

where  the  constant  depends  on  p  and  Cp+i  . 

Using  the  simple  estimates  for  H ?  in  terms  of  maximum  norms  gives 
Theorem  (2.4):  For  every  j  >  1  and  any  a  >  0,  we  can  choose  p  sufficiently  large 
so  that 

j— - rJ+l+O 

Hf{t)  <  const  •  (2.18) 

where  the  constant  depends  on  p  and  Cp+i  . 

Theorem  (2.1)  now  follows  from  Parseval’s  relation. 

It  may  seem  curious  that  the  initial  conditions  satisfy 

Hj  (0)  <  const  - - j55-, 

while  we  are  able  to  prove  that  (2.18)  holds.  However,  (2.18)  can  be  derived  from 
(2.17)  as  follows  (for  convenience  we  drop  the  a’s)  : 

Hj  <  const 
>  ~  vs 


=>  <  const  H*+i  <  const 


Z7u 


»+* 


r+J 


and  thus  using  the  interpolation  inequality  (2.13) 
l^’ulL  <  const 

ni.l*  )/(r- 1) 


<  const  (- 


<  const 


I/F+1 
•j  +  l  +  a 
oo _ 

v  j-l+a 


l)(|£>u|^)^"i)/^-l) 


\Dn\ 


(2.18)  now  follows. 


2.3  Estimates  for  two  space  dimensions 

In  this  section  we  obtain  a  sharp  bound  for  iDul^  for  two  dimensional  flows. 
In  this  case  the  incompressible  Navier-Stokes  equations  can  be  written  in  vorticity 
form 


6  +  +  v(y  =  ,  v  >  0,  (2.19a) 

+  vt  =  0  ,  ut  —  vm  =  £,  (2.19b) 

where  (  is  the  vorticity  and  u ,  v  are  the  velocity  components. 

Lemma  (2.2):  The  solutions  of  (£.19)  satisfy  the  maximum  principle 

l*(*,0loo  <  |*M)L. 

Proof:  This  well  known  result  follows  from  the  fact  that  at  a  local  maximum 
(minimum)  of  £*  =  =  0  and  <  0  (>  0). 

We  would  first  like  to  show  that  IZhif,*,  is  bounded  for  all  time.  Note  that  our 
energy  estimates  of  the  previous  section  are  still  valid  if  we  integrate  to  t  =  T  >  0 
instead  of  integrating  to  t  =  oo.  For  this  section  only,  let  us  redefine  the  quantities 
which  depend  on  this  bound  on  the  time.  For  example  we  define 

|0u|  =  supIDuU  and  Lr  =  f  Hl(t)dt. 

t<T  Jo 

We  know  from  basic  results  that  |Z?u|oo  exists  and  is  bounded  for  some  finite  time 
interval  [0,7’].  We  will  now  derive  estimates  for  jPuj^  which  are  independent  of 
T.  It  follows  from  the  results  of  the  previous  section  that  we  can  obtain  estimates 
for  all  derivatives  which  are  also  independent  of  T.  Then  from  well  known  results 
we  can  conclude  that  exists  and  satisfies  these  same  bounds  for  all  times. 

Lemma  (2.3):  For  any  aj  >  0  there  exists  a  constant  C(at)  such  that 

l*Hoo  =  BUP IDujoo  <  C(ai)|{(-,0)|io+o>i/-°1. 


(2.20) 


Proof:  For  any  /3  with  0  <  (3  <  1  we  define  the  Holder  semi-norm  by 


|^/|oo  =  8»P 


l/(*0  -  /(*»)! 

l*i  - 


Using  the  notation 

jD^uloo  =  max{]£>^u|oo,  |X>^w|oo} 

the  usual  Holder  estimates  for  the  solutions  of  Laplace’s  equation  (see  for  example 
[10]  ),  tell  us  that  for  any  (3  >  0  there  is  a  constant  C(/3)  such  that 

|D'-»uU  <  (2  21) 


Also,  the  convexity  of  Holder  norms  (see  [13]) 

\Dh+af\oo  <  con.t  |Dfcl+0l/|,00l^+°,/l»‘  , 

k  +  a  =  i(ki  +  Qi)  +  (1  —  <)(*!  +  02)  >  1,  0  <  i  <  1,  a,  ai ,  aj  >  0, 
and  Young’s  inequality  give  us  for  any  e  >  0 

|Z?uU  <  e|Z?lu|00  -f  const  c~0\Dl~0u\oc- 


Using  the  Sobolev  inequality  (2.12)  for  |i?2u|0O  and  (2.21)  we  obtain 

tT/s 


\D*L  <  conti  e^jf  -  +  const  (~p 


t7/’ 


<  con.t  [e!^k-+e^C(^|e 


Choosing 


<■+«  =  cw\i\ 


_  „'/! 


00 - 7/ J  ' 

\DU\1 


\Du\„<  con.t  C(/S)|(L  fc(g)KL=^7?  1 


0 

TTB 


gives 


Thus 


i^uioc  ^  c°ntt  (cmtL) 


—  xl/(l-60/») 


-70/(t-S0) 


and  the  lemma  follows  since  <  |^(-I0)Joo. 

In  two  space  dimensions  estimates  on  the  vorticity  appear  more  naturally.  In 
[17]  we  proved  the  results  of  theorem  (2.1)  in  the  two-dimensional  case  using  the 
vorticity  formulation  of  the  equations.  In  that  paper  the  quantities 


«0"’  +  Ol’ 

take  the  place  of  the  H *.  The  estimate  corresponding  to  (2.18)  is 


J*(t)  <  const  ffp+\(t)  <  const  — 


(2.22) 


(2.23) 


88 


v»>y/ 

» 


fijfe 


We  refer  to  the  J„  in  the  section  on  numerical  results. 


I 


3  Numerical  Results 

We  first  describe  the  procedure  we  use  to  numerically  solve  the  two-dimensional 
Navier-Stokes  equations.  In  brief,  we  discretize  in  space  using  the  Fourier  (pseudo- 
spectral)  method  and  solve  in  time  using  a  fourth  order  predictor-corrector  method. 
The  equations  are  solved  in  Fourier  space  and  the  diffusion  term  i/Au  is  treated  in 
a  fully  implicit  manner.  We  now  proceed  to  present  more  details. 

We  solve  the  two-dimensional  incompressible  Navier-Stokes  equations  in  the 
vorticity  stream  function  formulation: 


it  +  («0»  +  (r0v  =  +  / 

A1>  = -i,  («,»)  =  (V’m-V’a)- 


(3.1a) 

(3.1b) 


The  computational  domain  is  taken  to  be  a  2k  periodic  square.  The  solution  is 
represented  as  a  truncated  Fourier  series  with  u  denoting  the  discrete  approximation 
to  £  and  u)  denoting  the  discrete  Fourier  transform  of  u>: 

»(«,»,«)=  Y.  E  *(*••*». 

Similarly  the  Fourier  transform  of  and  /  are  denoted  by  i>  and  /  respectively. 
The  equation  for  the  Fourier  coefficient  u>(fci ,k2,t)  is 


ut  +  iki(uv)  4-  ik2(vu>)  =  -i/(fcf  +  fcj)u>  4-  / 
(k{  4-  kl)j,  =  & 


(3.2a) 

(3.2b) 


The  convolutions  uu>  and  vu>  (i.e.  the  Fourier  transforms  of  the  products  uu>  and  vw) 
are  computed  from  a,  v  and  u>  by  transforming  to  real  space,  forming  the  products 
and  then  transforming  back  to  Fourier  space,  (pseudo-spectral  method).  It  is  not 
hard  to  see  that  the  computation  of  o>,  can  be  done  with  five  two-dimensional  fast 
Fourier  transforms  (FFT’s).  In  fart,  only  four  FFT’s  are  needed  since  one  can  write 
(reference  Basdevant  [2]) 

=  ((iUa  -  -  (V^,),,)- 


However,  in  the  calculations  presented  here  the  less  efficient  method  was  used. 

The  equations  (3.2)  can  be  written  in  the  form  of  a  large  system  of  ordinary 
differential  equations: 

^  =  F(y,()  (3.3) 

where  y  is  the  vector  with  components  u>(ki,ki). 

Time  stepping  is  performed  using  a  predictor-corrector  applied  directly  to  equa¬ 
tion  (3.3).  Let  yn  denote  the  approximation  to  y(nA<)  and  Fn  =  F(yn,nAt).  We 
use  the  fourth  order  Adams  predictor-corrector  scheme  given  by 

At 

y,  =  yn  +  y^-(23Fn  -  16Fn_!  +  5Fn_2)  (3.4a) 

y»  +  l  =  yn  +  ^(9FP  +  19Fn  -  5Fn_i  +  Fn_2).  (3.4b) 

Here  yp  is  the  result  of  the  Adams-Bashforth  predictor,  Fi>  =  F(yj>.(»  +  !)A0 
and  yn+i  is  the  corrected  value  obtained  from  approximating  the  implicit  Adams- 
Moulton  scheme.  A  single  time  step  thus  requires  two  evaluations  of  the  right  hand 
side  F.  The  classical  fourth  order  Runge-Kutta  method  is  used  to  obtain  starting 
values  for  (3.4).  These  are  required  initially  and  whenever  the  time  step  is  changed. 

For  stability  reasons  one  may  want  to  integrate  the  diffusion  term,  uAu,  in  an 
implicit  manner.  In  the  Fourier  representation  this  term  is  very  simple  and  thus  can 
be  easily  treated  in  a  fully  implicit  and  accurate  manner.  We  write  the  equations 
(3.3)  in  the  following  way 

-Jl  =  G(y.O  -  Ay  A  =  diag(..., v(k\  +  fe2),...). 

where  the  right  hand  side  F  has  been  split  with  A  the  diagonal  matrix  corresponding 
to  the  diffusion  term.  This  last  equation  can  be  written  in  the  form 

|(*AV)  =  eA'G. 

Now  apply  the  time  stepping  procedure  (3.4)  to  this  equation  viewed  in  terms  of  the 
new  dependent  variable  eA<y.  After  division  by  eA<  the  predictor-corrector  scheme 
which  results  is 

y p  =  e-AA<y„  +  ^(23e-AA‘Gn  -  lee'^'G,^  +  5e"SAA,Gn_2)  (3.5a) 
y«  +  i  =  e-A*‘y»  +  ^(9G,  +  19e-AA,Gn  -  5e-2AatGn_!  +  e-SA*‘Gn_2).(3.5b) 


The  Runge-Kutta  scheme  is  transformed  in  a  similar  fashion.  The  terms  e~^k>  +kl)At 
are  stored  and  need  only  be  recalculated  when  At  changes.  These  resulting  schemes 
are  exact  in  the  absence  of  the  convection  terms  (G  =  0). 

The  variable  time  step  is  chosen  by  stability  and  accuracy  considerations  with 
At  chosen  to  satisfy  the  condition 

CFLmi„  <  (Noo  +  Moo)^  <  CFLmmJt  (3.6) 

h 

where  h  =  2 /IV,  ( N  =  max(!Vi ,  JVj)).  The  stability  region  of  the  explicit  predictor- 
corrector  method  (3.4)  is  shown  in  figure  1  . 


Figure  1  Stability  region  for  the  predictor-corrector  scheme. 

When  (3.4)  is  applied  to  the  model  problem  y'  =  Ay,  the  time  step  is  restricted 
by  (approximately)  |A|Af  <  1.2  if  A  is  purely  imaginary  and  by  — AAt  <  1.9  if  A  is 
real.  One  expects  the  implicit  predictor-corrector  scheme  (3.5)  to  have  better  sta¬ 
bility  properties  than  the  explicit  one  (3.4).  CFLm;n  and  CFLmkX  are  the  minimum 
and  maximum  allowable  values  for  the  Courant-Friedrichs-Lewy  number: 


CFLmax  would  be  taken  less  than  the  stability  limit  for  the  model  problem.  (The 
choice  of  h  in  our  definition  of  CFL  instead  of  the  true  h  =  2ir/N  means  that  we  can 
compare  CFL  directly  to  the  normal  stability  limit  for  the  model  problem.)  When 
the  condition  (3.6)  is  violated  the  new  time  step  is  chosen  so  that 

(Moo  +  Moo)^  =  CFLopt. 
h 

3.1  Verification  of  the  Numerical  Approximation 

In  this  section  we  present  results  which  illustrate  the  accuracy  of  the  numerical 
approximation  that  we  use.  In  test  1  we  show  that  the  time  stepping  procedure 
is  accurate  to  fourth  order  in  At.  In  test  2  we  consider  the  convergence  of  the 
numerical  solution  as  the  number  of  modes  is  increased. 

Test  1:  Accuracy  of  the  time  stepping  procedure 

It  is  easy  enough  to  choose  the  forcing  /  in  the  Navier-Stokes  equations  (3.1) 
so  that  the  true  solution  is  known  to  be  some  given  function.  Numerous  tests  of  this 
kind  were  performed.  In  all  cases  the  numerical  solutions  converged  to  the  exact 
solutions  at  a  rate  very  close  to  fourth  order  in  the  time  step  At. 

As  a  more  realistic  study  of  the  convergence  of  the  time  stepping  routine  we 
consider  a  sequence  of  calculations  with  fixed  random  initial  data  and  decreasing 
time  steps.  The  initial  conditions  are  identical  with  those  used  in  section  (3.2)  for 
the  decay  of  random  initial  data.  Keeping  the  same  initial  conditions,  and  with 
N\  =  Nt  =  128,  v  =  10-4,  the  equations  were  solved  with  three  different  (fixed) 
time  steps:  At  =  .05,  At  =  .025  and  At  =  .0125.  The  computed  maximum  value 
for  the  CFL  number  in  each  run  was  1.2,  .6  and  .3,  respectively.  (Recall  that  the 
stability  limit  for  the  explicit  version  of  the  predictor-corrector  scheme  is  about  1.2 
on  the  imaginary  axis.  We  were  able  to  obtain  good  results  for  values  of  the  CFL 
number  as  large  as  1.5,  which  substantiates  the  belief  that  the  implicit  predictor- 
corrector  scheme  has  better  stability  properties.)  We  use  the  results  from  the  three 
runs  to  estimate  the  rate  of  convergence  as  a  function  of  At  as  well  as  to  estimate 


the  actual  error.  We  measure  both  the  discrete  maximum  error  and  the  l2  error 
defined  by 


M«  =  max and  \u\\  :=  — —  V  |w(*i,yj)|2. 

(»i.Vy)  1  ■**  2  .  . 


(*(.»>) 


By  assuming  that  the  computed  solution  is  converging  to  the  true  solution  as  0(Atp) 
we  can  determine  approximate  values  for  p  and  the  errors 


£«,(<,  At)  :=  |u/computed(-,t;  At)  -  «tru ,(•,<;  Af)U  =  0(Atp), 
£2(t,  At)  :=  |u>computed(-,  t;  At)  -  wltue(-,f;  Af)|,  =  0(Afp). 


These  values  are  given  in  table  I  for  the  maximum  norm  errors  and  in  table  II  for 
the  l2  errors. 


E00(At  =  .05) 

Eoo(At  =  .025) 

^oo(At  =  .0125) 

0.68  x  10“2 

0.47  x  10"8 

0.32  x  10-4 

0.27  x  10~2 

0.18  x  10-8 

0.13  x  10~4 

0.97  x  10“s 

0.69  x  10“4 

0.49  x  10-6 

0.69  x  10“s 

0.47  x  10~4 

0.32  x  10~s 

0.48  x  IQ"3 

0.31  x  10-4 

0.20  x  10-s 

Table  II  -  Estimated  maximum  errors  and  convergence  rate:  0(Atp) 


t 

#2  (At  =  .05) 

^2  (At  =  .025) 

E2{At  =  .0125) 

0.40  x  10 
0.17  x  10 
0.73  x  10 
0.56  x  10 
0.48  x  10 


0.29  x  10~4 
0.12  x  10~4 
0.50  x  10“5 
0.37  x  10-5 
0.31  x  10”5 


0.21  x  10-* 


0.86  x  10 


0.34  x  10 


0.24  x  10 


0.20  x  10 


Table  II  -  Estimated  l2  errors  and  convergence  rate:  0(Atp) 


Test  2:  Convergence  for  random  initial  data 

We  consider  the  computation  which  is  described  in  section  (3.2)  under  the 
heading  of  Run  1:  Decay  of  Random  Initial  Data.  This  computation  was  run  with 
N  =  N\  =  N2  =  128  (u>i28)  and  also  with  N  =  256  (u»2ss)-  The  initial  conditions  for 


vv 

A- 

<V 


mm 


the  two  runs  were  the  same  to  single  precision  (about  6  —  7  decimal  digits),  although 
the  actual  computations  were  done  in  double  precision.  The  variable  time  step  was 
determined  by  the  parameters  (CFLmin, CFLoptl CFLmaJt)  =  (.8,1.,  1.2)  In  table  II 
we  indicate  the  maximum  difference  and  the  Ij  difference  between  the  two  runs  at 
various  times.  Due  to  the  variable  time  step  the  solutions  were  not  compared  at 
exactly  the  same  times.  The  difference  between  the  times  is  given  in  the  table  as 
tj66  —  tus-  Note  that  the  maximum  difference  between  the  two  solutions  occurs  at 
smaller  times.  Later  on  when  the  solution  becomes  smoother  the  errors  are  smaller. 
Further  details  of  this  run  can  be  found  in  the  next  section. 


l^ase  oo 

|w268  “  «128|oo 

|w2se  -  wus|2 

<2SS  —  <128 

1.00 

.10  x  10-s 

.38  x  10~8 

0.0 

.75 

.55  x  10-1 

.45  x  10-2 

+.42  x  10"5 

.67 

.20  x  10"1 

.29  x  10"2 

-.21  x  10~8 

.60 

.15  x  10"1 

.21  x  10~2 

+.23  x  10"3 

.59 

.10  x  10"1 

.73  x  10~8 

+  .18  x  10“4 

.58 

.73  x  10“J 

.58  x  10"* 

+.23  x  10-4 

.57 

.43  x  10-2 

.88  x  10-8 

-.19  x  10-8 

.56 

.51  x  10"2 

.88  x  10~8 

+.25  x  10-8 

.56 

.26  x  10"2 

.38  x  10"8 

+  .37  x  10-4 

.56 

.27x10-* 

.36  x  10~8 

+  .12  x  lO"4 

.55 

.31  x  IQ’2 

.37  x  10"8 

-.62  x  lO"4 

'Ih.ble  II  -  Convergence  of  random  initial  data  v  =  10 


3.2  Computational  Results 

In  this  section  we  present  the  results  of  four  different  runs: 

(1)  Run  I:  Decay  of  random  initial  data,  v  =  10-4,  N  =  256  and  N  =  128. 

(2)  Run  II:  Decay  of  random  initial  data,  v  =  10“ *,  N  =  512. 

(3)  Run  III:  Decay  of  smooth  random  initial  data,  i/=2x  10”*,  N  =  256. 

(4)  Run  IV:  Random  forcing,  i/  =  .5x  10-s,  v  —  .5  x  10“4. 

Run  I:  Decay  of  random  initial  data,  v  —  1 0 ~ 4 ,  N=256  and  N=128. 

For  the  first  run  we  consider  the  time  evolution  of  the  Navier-Stokes  equations 
for  random  initial  data.  The  initial  conditions  for  the  vorticity  were  chosen  60  that 

=  (k  +  ^jtyy  *  = 

with  a  random  phase.  (Actually  the  initial  spectrum  was  set  to  sero  for  all  wave 
numbers  above  some  large  value  of  k.)  The  constant  C  was  determined  by  nor¬ 
malizing  the  maximum  value  of  the  vorticity  to  be  1  at  t  =  0,  )«;(•,  *,  O)!*,  =  1. 
The  value  of  the  viscosity  was  taken  as  v  =  10-4  and  the  number  of  modes  was 
—  N2  =  256.  We  show 

1)  contour  plots  of  the  vorticity  (figure  2).  Dashed  lines  indicate  negative  contour 
values. 

2)  surface  plots  of  ,  kj)  in  the  cosine-sine  representation.  The  discrete  Fourier 
series  for  u>  is  actually  represented  in  the  computer  code  as  a  real  series  in 
cosines  and  sines.  The  surface  plot  shows  the  magnitude  of  the  coefficients  of 
this  series.  The  coefficients  are  ordered  in  the  following  manner: 


CiCi 

Cl«l 

C1C2 

Cl  *2 

ClC8 

«lCl 

*1*1 

*lC2 
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CjCl 

C2«1 
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C2C3 
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*2*1 

*2®2 

*2*2 

*2Cs 

C3C1 

Cj*i 

C3C2 

CS»2 

Cs  C3 

where  c*cj  is  the  coefficient  of  cos(kz)  cos(ly),  c*sj  the  coefficient  of  cos(fcz)  sin(ly) 
and  so  on.  The  lowest  frequency  modes  are  located  at  the  top  of  the  surface 
plot  (figure  3  ).  Only  the  first  128  modes  are  shown  in  the  surface  plots. 


3)  plots  of  the  energy,  enstrophy,  Jj,  Jj,  as  a  function  of  time,  and  the  decay 
of  the  vorticity  spectrum  as  a  function  of  k.  In  figure  4  a  we  plot  the  square 


root  of  the  total  energy 


the  square  root  of  the  enstrophy 


HI. 


and 

=  w'/’dKII1  +  IKII1)1'1 

as  functions  of  time.  In  figure  4  b  we  plot  the  normalized  versions  of  Jx(<), 
i/2/2  Jj(<)  and  Recall  that 

j;<(> = O’ + O'- 


In  each  case  the  functions  plotted  are  scaled  so  their  maximum  value  is  1.  This 
maximum  value  is  indicated  on  the  plot  as  the  value  of  Scale.  In  figure  4  c  some 
selected  Fourier  coefficients  are  plotted  as  functions  of  time.  Finally  in  figure 
4  d  we  show  log-log  plots  of  u(k)  versus  k.  The  quantity  u>(k),k  =  1,2,... 
is  defined  to  be  the  average  value  of  |<I>(/i,/2)|  over  all  wave  vectors  (f x , )  for 
which  k  is  the  closest  integer  to  l  =  |(/i,/j)|: 

*(*)=(  £  M<..Ml)/(  £  i 

\i-ai<i/i  '  V-* i<i/i 

We  plot  logi0(w(fc))  versus  log10(k)  for  different  times.  Lines  with  slopes  -1 
and  —2  are  also  marked.  Note  that  if  ~  k~a  then  E{k)  ~  fc_2a_1. 


For  comparison,  in  figures  5  -  6  ,  we  show  the  results  of  the  same  run  when  only 
half  as  many  modes  were  used,  N  =  128.  Essentially  the  only  noticeable  difference 
is  in  the  plot  of  the  spectral  decay.  A  quantitative  comparison  of  the  256  and  128 
runs  was  given  in  section  (3.1). 
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Run  II:  Decay  of  random  initial  data,  iv=10“s,  N  =  512. 
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The  value  of  v  is  taken  as  10-s.  The  initial  conditions  are  the  same  as  Run  I. 
The  number  of  modes  was  JVj  =  Nj  =  512.  The  results  are  shown  in  figures  7  - 
8  .  Note  that  for  technical  reasons  the  contour  plots  were  made  by  projecting  the 
solution  to  a  256  x  256  grid.  It  is  interesting  to  compare  this  run  (i/  =  I0~s)  to  the 
previons  run  (v  =  10~ "*). 

Run  III:  Decay  of  smooth  random  initial  data,  v  —  2  x  10~s. 

In  this  run  we  begin  with  initial  data  which  is  much  smoother  than  in  the 
previous  runs.  The  initial  vorticity  spectrum  is  chosen  so  that 

\u>(kuk2)\  =  Ck*-tth/h°f  ,  fco  =  3.5 

with  random  phase.  The  constant  C  is  chosen  so  that  |w(-,  •,  0)1^  =  1.  The  viscosity 
was  2  x  10-s  and  the  number  of  modes  was  N  =  256.  These  initial  conditions  are 
similar  to  those  used  by  Brachet  and  Sulem  [7].  We  have  run  for  longer  times  than 
the  results  shown  in  [7j.  Plots  for  this  run  are  given  in  figures  9  -  10  . 

Run  IV:  Random  forcing. 

In  this  run  we  consider  the  problem  when  the  equations  are  forced  in  a  range  of 
low  Fourier  modes.  For  the  forced  problem  there  appears  to  be  no  easy  way  to  obtain 
a  sharp  bound  on  the  maximum  of  the  vorticity.  We  have  found  experimentally  that 
when  the  forcing  is  chosen  to  be  0(1)  the  solution  grows  and  does  not  remain  0(1). 
For  example  in  figure  11  we  show  the  results  of  a  run  in  which  the  forcing  /  is 
chosen  bo  that  |/|oo  =  1  and  in  which  the  initial  vorticity  is  rero.  In  particular  the 
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where  the  scaling  factor  C,  was  chosen  to  ensure  that  |/|oo  =  1-  In  this  run 
u  =  10~8,  Ni  =  N3  -  256  and  (CFLmin,  CFLopt,  CFLm«x)  =  (.5,  .8,1.)-  In  figure 
11  b  we  have  made  plots  of  the  maximum  norms 

MAX(U)  =  max(|tt(-,,t)|OO)|»(I10loo) 

MAX(DU)  =  max(|^(-,  •,*)!<»,  I^(‘,  ‘•OU.  I^(‘.  *,0U,  I^O,  ’.OU). 

and 

MAX(W)  =  |«(-,-,l)U. 

Even  by  time  t  —  200  the  solution  continues  to  grow. 

In  contrast  when  the  forcing  is  chosen  to  be  0( i/)  we  do  not  see  growth  in  the 
)<*>(•,  •,/)|00.  This  observation  is  presented  in  figures  12-13  where  we  have  made 
runs  with  u  =  .5  x  10~8  ( N  =  128)  and  v  =  .5  x  10-4  ( N  =  256).  The  initial 
conditions  for  |u/|  were  defined  by  the  matrix  of  coefficients  given  above  but  in  this 
case  the  constant  C,  was  chosen  so  that  |u>(-,  •,  0)joo  =  1.  The  forcing  was  constant 
in  time  and  defined  from  the  relation 

i/Au  +  f  ~  0. 


4  Discussion 

We  have  shown  that  for  both  two  and  three  dimensional  flows  the  minimum 
scale,  Amjn,  is  essentially  proportional  to  i/1^,/|Du|oc  ^  .  We  now  relate  this 
minimum  scale  to  the  decay  rate  of  the  energy  spectrum. 

Let  us  assume  that  at  a  given  time  t  the  energy  spectrum  has  a  power  law 
behavior  in  some  range  of  wave  numbers,  the  inertial  range, 

E(k)  ~  k !  =  0(1)  <  k  <  0( l/\min). 

We  have  proved  that  the  quantities 


remain  of  order  one,  provided  they  are  initially  so.  When  this  order  one  bound  is 
achieved  we  call  the  flow  maxima/  dissipative.  Note  that  2vH\  is  the  rate  of  energy 
dissipation,  e,  and  that  the  rate  of  enstrophy  dissipation  rj  is  bounded  by  vH\. 

Assuming  that  the  leading  order  contribution  to  the  integral  for  Hj(t)  is  deter¬ 
mined  from  wave  numbers  in  the  inertial  range  and  using  the  power  law  behavior 
for  E(k)  we  obtain 


i/p-x  ,  i/p_1  .  I  Du 

— ~ +T*,(0  ~  =-7TT  /  klpE(k)dk  ~  I - 

\Du\’+  |Dul  +1  Jhi  v 


-0-1/2 
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■0  +  8/1 


In  two  space  dimensions  we  know  that  l-Dul^  is  essentially  bounded  by  the 
maximum  norm  of  the  initial  vorticity.  Let  us  thus  assume  that  the  initial  values 
are  scaled  so  that  l-Dul^  is  of  order  one.  In  this  case 

=^Tr",’< I) (4.1) 

|0»L 

and  for  maximal  dissipative  flows  it  follows  that  E(k)  ~  fc~8,  the  power  law  behavior 
predicted  by  the  Batchelor-Kraichnan  theory  [3]  [16]. 

In  three  dimensions  if  we  speculate  that  jZTuj^  =  0{v~1)i  then 


_ 17*  ^  7(0  +  l/*)+0-S/l 

F  +  l^P 


In  this  way  for  maximal  dissipative  flows,  we  obtain  a  relation  between  the  power 
law  behaviour  of  the  energy  spectrum  and  the  sise  of  IZhil^: 

2  +  2  7 

When  l-Dul^  =  0( and  7  =  1/2,  we  obtain  0  =  5/6  and  E(k)  =  &-5/8,  the 
power  law  behavior  predicted  by  Kolmogoroff  [15]. 

We  now  return  to  the  numerical  results  of  the  random  initial  data  runs,  refer¬ 
ence  figures  2-4  and  7  -  8  .  The  initial  conditions  were  chosen  so  that  E[k)  ~  k~s. 
The  numerical  results  show  that  this  k~3  power  law  seems  to  remain  over  an  initial 
time  interval.  The  numerical  results  further  indicate  that  as  the  flow  evolves,  the 


quantities  upJ which  behave  like  v*Hp+i,  slowly  decrase  and  the  energy  spectrum 
steepens.  In  a  later  regime  the  flow  is  dominated  by  the  presence  of  large  regions 
of  relatively  constant  vorticity.  There  seems  to  be  some  evidence  from  the  figures 
to  suggest  that  the  quantities  J*  decay  by  a  factor  of  about  Amjn  =  and  that 
t(k)  ~  *-ls,  E{k )  ~  *"4.  When 
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the  argument  in  (4.1)  predicts  f3  =  2,  and  is  thus  consistent  with  the  numerical 
results.  These  results  are  in  agreement  with  Saffman’s  theory  for  two  dimensional 
turbulence  [22]. 

In  figure  1  we  outline  an  hypothesized  behaviour  of  vvJ*(t)  for  the  decay  of 
two-dimensional  turbulence.  In  the  first  stage  of  development  of  the  flow,  vv  J *  may 
show  an  overall  increase  as  the  flow  evolves  to  a  state  of  maximal  dissipation.  (Of 
course,  depending  on  the  initial  data,  this  maximal  dissipative  state  may  never  be 
reached.)  This  dissipation  rate  can  not  continue  for  a  long  time  interval  but  must 
decrease.  The  power  law  then  slowly  changes  from  fc~8  to  a  more  rapid  decay.  The 
flow  becomes  organized  into  coherent  structures,  a  regime  with  iSJ*  ~  vx I7  (?) 
and  where  Saffman’s  theory  would  predict  E(k)  ~  fc-4.  This  regime  presumably 
exists  for  long  times,  since  the  viscosity  now  plays  a  minimal  role.  This  scenario  is 
suggested  by  our  computations  and  other  similar  ones.  In  particular,  Brachet  and 
Sulem  [7]  show  high  resolution  computation  with  initial  data 

E(k)  ~  cke~(k/ho'>\ 

similar  to  the  one  presented  in  this  paper  (reference  figures  9  -  10  ).  They  found  an 
increase  in  the  energy  power  law  reaching  a  maximum  at  about  fc-s.  At  this  stage 
the  rate  of  enstrophy  dissipation  is  maximum,  in  accordance  with  our  analytical 
results. 


It  is  conceivable  that  a  similar  scenario  is  present  in  the  decay  of  three  dimen¬ 
sional  turbulence.  In  two  dimensions  when  large  coherent  structures  are  formed 
the  main  contribution  to  Hp{t)  comes  from  one  dimensional  layers  separating  these 


width  Amjn.  Thus  in  either  case  we  obtain  a  new  relation  between  the  power  law 
behaviour  of  the  energy  spectrum  and  the  size  of  l-Du^  =  0(i/-7): 
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Figure  13  Contour  plots  of  vorticity  spectrum  for  random  forcing.  '  f 
Vr.5v  10-3,  N  =  128  and  v  -  .5  *  1 0 "  ^ ,  X  -  25C 
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